---
title: "R8.SoilType"
output: pdf_document
---

# Load Packages
```{r}
library(tidyverse)
library(stargazer)
library(gridExtra)
```

# Load Data
```{r}
load("Data/points_soiltype.Rdata")
```

# Run Regressions Comparing Soil Type A to Soil Types B, C, D, and E (App. Table 3aii)
```{r}
m.1 <- lm(log(Av_Pt_Ha_2) ~ log(Av_Pt_Ha_1), points_soiltype)
m.2 <- lm(log(Av_Pt_Ha_3) ~ log(Av_Pt_Ha_1), points_soiltype)
m.3 <- lm(log(Av_Pt_Ha_4) ~ log(Av_Pt_Ha_1), points_soiltype)
m.4 <- lm(log(Av_Pt_Ha_5) ~ log(Av_Pt_Ha_1), points_soiltype)

stargazer::stargazer(list(m.1, m.2, m.3, m.4), type = "text", 
                     covariate.labels = c("Log(Av. Soil A)"), 
                     dep.var.labels = c("Log(Av. Soil B)", "Log(Av. Soil C)", "Log(Av. Soil D)", "Log(Av. Soil E)"))
```

# Plot Relationship Between Soil Types (App. Table 3aiii)
```{r}
g.1 <- ggplot(points_soiltype, aes(x=log(Av_Pt_Ha_1), y=log(Av_Pt_Ha_2))) + 
  geom_point() + geom_smooth(method='lm', formula = y ~ x) +
  theme_light() + xlab("Log(Average Points/Ha, Soil Class A)") + ylab("Log(Average Points/Ha, Soil Class B)")
g.2 <- ggplot(points_soiltype, aes(x=log(Av_Pt_Ha_1), y=log(Av_Pt_Ha_3))) + 
  geom_point() + geom_smooth(method='lm', formula = y ~ x) +
  theme_light() + xlab("Log(Average Points/Ha, Soil Class A)") + ylab("Log(Average Points/Ha, Soil Class C)")
g.3 <- ggplot(points_soiltype, aes(x=log(Av_Pt_Ha_1), y=log(Av_Pt_Ha_4))) + 
  geom_point() + geom_smooth(method='lm', formula = y ~ x) +
  theme_light() + xlab("Log(Average Points/Ha, Soil Class A)") + ylab("Log(Average Points/Ha, Soil Class D)")
g.4 <- ggplot(points_soiltype, aes(x=log(Av_Pt_Ha_1), y=log(Av_Pt_Ha_5))) + 
  geom_point() + geom_smooth(method='lm', formula = y ~ x) +
  theme_light() + xlab("Log(Average Points/Ha, Soil Class A)") + ylab("Log(Average Points/Ha, Soil Class E)")

g.all<-grid.arrange(g.1, g.2, g.3, g.4, ncol = 2)

ggsave(plot = g.all, "Figures/Afg3a_PointClass.pdf", width = 8, height = 8)
```